************Here is the data set up for
the Nonlinear SEM example 2 for the interaction model in Wall (2009);

   data param;

    reliab = .75;
	rsquare_seminter = .5;
              gam0= .2 ; gam1= .1 ; gam2= .1 ; gam3= .2;
              b10 = 0 ; b20 = 0 ; b30 = 0 ; b40 = 0 ; b50 = 0 ; b60 = 0 ;
              b11 = .7 ; b21 = .3 ;
              b32 = .4 ; b42 = .8 ;
              b53 = .5 ; b63 = .4 ;

              mf1 = 5; mf2 = 2;
              ph1= 1 ; ph12=.5;
                       ph2 = 1;
              ****zeta= 11 ;  **** this zeta is defined below to govern the R^2 of the interaction equation which I've done with the rsquare_seminter;

	   varf1f2 = (ph1+mf1*mf1)*(ph2+mf2*mf2) + (ph12+mf1*mf2)**2 - 2*mf1*mf1*mf2*mf2;
       cov1f1f2= mf1*ph12+ph1*mf2;
       cov2f1f2= mf2*ph12+ph2*mf1;

	   noeqerror = gam1*gam1*ph1 + gam2*gam2*ph2 + gam3*gam3*varf1f2
               + 2*gam1*gam3*cov1f1f2 + 2*gam2*gam3*cov2f1f2 + 2*ph12; 
       zeta = noeqerror*(1-rsquare_seminter)/rsquare_seminter;

       phf3 = noeqerror + zeta;
	   *phf3 =7.9; ***this was the variance of f3 when the data was uniform, but when it is normal the formulas work;


       ps1 = (b11*b11*ph1*(1-reliab)/reliab);
       ps2 = (b21*b21*ph1*(1-reliab)/reliab);
       ps3 = (b32*b32*ph2*(1-reliab)/reliab);
       ps4 = (b42*b42*ph2*(1-reliab)/reliab);
       ps7 = (ph1*(1-reliab)/reliab);
       ps8 = (ph2*(1-reliab)/reliab);

       ps5 = (b53*b53*phf3*(1-reliab)/reliab);
       ps6 = (b63*b63*phf3*(1-reliab)/reliab);
       ps9 = (phf3*(1-reliab)/reliab);

meanoff3 = gam0 + gam1*mf1 + gam2*mf2 + gam3*(mf1*mf2+ph12);

     run;
proc print data = param; run;

        proc iml;
        use param;
        read all var{ph1 ph2 ph12} into temp;
        phi=j(2,2,0); out=j(1,3,0);
        phi[1,1]=temp[1,1];
        phi[2,2]=temp[1,2];
        phi[1,2]=temp[1,3];
        phi[2,1]=temp[1,3];

        halfphi= half(phi);
        out[1,1]= halfphi[1,1];out[1,2]= halfphi[1,2]; out[1,3]= halfphi[2,2];
       create param1 from out[colname = {hph11 hph12 hph22}];
       append from out;
       close param1;
       quit;
     data param2; merge param param1; run;
proc print data = param2; run;
proc datasets nolist;
     delete  saveparm save;
run;

%macro dosimulate(sam,obs);
     data start;
      do i=1 to %eval(&sam*&obs);
         xf1=rannor(125);  /*xf1=ranuni(125); xf1=sqrt(12)*(xf1-.5)*/;
         xf2=rannor(125);  /*xf2=ranuni(125); xf2=sqrt(12)*(xf2-.5)*/;
         %do j=3 %to 12;
         x&j=rannor(125);
         %end;
         output;
      end;

     data a;  set start;
     if _n_=1 then set param2;

     f1= mf1 + hph11*xf1;
     f2= mf2 + hph12*xf1 + hph22*xf2;
     f3= gam0 + gam1*f1 + gam2*f2 + gam3*f1*f2 + sqrt(zeta)*x3;
     Y1   = b10+ b11*f1   + sqrt(ps1)*x4;
     Y2   = b20+ b21*f1   + sqrt(ps2)*x5;
     Y3   = b30+ b32*f2   + sqrt(ps3)*x6;
     Y4   = b40+ b42*f2   + sqrt(ps4)*x7;
     Y5   = b50+ b53*f3   + sqrt(ps5)*x8;
     Y6   = b60+ b63*f3   + sqrt(ps6)*x9;
     Y7   =          f1   + sqrt(ps7)*x10;
     Y8   =          f2   + sqrt(ps8)*x11;
     Y9   =          f3   + sqrt(ps9)*x12;
output;
     keep Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 Y9 f1 f2 f3;
     run;

     data a; set a;
          %do p=1 %to &sam;
          if %eval(1+(&p-1)*&obs)<=_n_<=%eval(&p*&obs) then iter=&p;
		  count = _n_;
		  dummy = 1;
          %end;
		  run;

          proc means data=a noprint; by iter;
          var y1-y9;
          output out=outmean mean=mean1-mean9;
          run;

          data a; set a;
           merge a outmean; by iter;
          run;

%mend dosimulate;

%dosimulate(1,500);

/*data b; set a (where = (iter = 1));
file "C:\file\dataexpb.dat";
put dummy count y1-y9 iter;
run;*/

/*
*Here I am putting just the y1-y9 since Winbugs needs to only have the input data
data and not other stuff;
data b; set a (where = (iter = 1));
file "C:\file\dataexpb_forwinbugs.dat";
put y1-y9;
run;*/